Two-Dimensional Diffusion in the Presence of Topological Disorder 
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How topological defects affect the dynamics of particles hopping between lattice sites of a dis- 
torted, two-dimensional crystal is addressed. Perturbation theory and numerical simulations show 
that weak, short-ranged topological disorder leads to a finite reduction of the diffusion coefficient. 
Renormalization group theory and numerical simulations suggest that longer-ranged disorder, such 
as that from randomly placed dislocations or random disclinations with no net disclinicity, leads to 
subdiffusion at long times. 

PACS numbers: 0.5.40Jc, 61.72.Lk, 66.30-h 



I. INTRODUCTION 

Diffusion in random media is a well-studied problem 
The mean-square displacement of a tracer particle 
behaves at long times in a way that depends on the char- 
acter of the random forces induced on the tracer by the 
disorder. Forces that arise from random potentials lead 
to a reduction of the transport, with subdiffusion pos- 
sible for diffusion of an ion in a medium with quenched 
charges obeying bulk charge neutrality 1] . Interestingly, 
the same subdiffusion results from diffusion of an ion 
in a medium with randomly-placed, quenched dipoles 
07 0,0, 0| . Forces that arise from entrainment along fluid 
streamlines lead to an increase in the transport, with the 
well-known result of turbulent super-diffusion possible for 
random streamlines with statistics characteristic of fluid 
turbulence 

Distortion of the underlying lattice upon which the dif- 
fusion occurs is a very different type of disorder. In par- 
ticular, topological defects such as dislocations or discli- 
nations should affect the transport properties of a dif- 
fusing tracer particle. These topological defects cause 
a global rearrangement of the connectivity of the lat- 
tice upon which the diffusion occurs. Moreover, there is 
an elastic response of the lattice to such defects, and so 
there is also local expansion or compression of the crystal 
unit cells. Study of how such topological defects affect 
the transport is, therefore, an interesting and challenging 
problem. Among other results, it might be expected that 
randomly-placed dislocations and random disclinations 
with no bulk disclinicity will lead to similar dynamics, 
given the results regarding dynamics in random poten- 
tials and the analogy between linear elasticity theory and 
electrostatics. 

Previous work has begun to address the question of 
how topological disorder affects the transport. Random 
disclinations, with no net disclinicity, were predicted to 
lead to subdiffusion [|| . A single dislocation, on the other 
hand, was predicted to increase the local diffusivity 0. 
These studies, however, were approximate 8J. In partic- 
ular, rotational symmetry was assumed in the dislocation 
problem, and no effects of lattice expansion or contrac- 
tion were allowed in the disclination problem. 



Transport in a two-dimensional crystal with topolog- 
ical defects, then, remains an interesting and unsolved 
problem. Our model of surface diffusion, and the Fokker- 
Planck equation that results, is introduced in Sec. II. 
How the topological defects affect the transport, and a 
field theoretic description used to analyze the dynamics, 
is described in Sec. III. Perturbation theory and com- 
puter simulation are used to examine the effect of non- 
singular topological disorder on the diffusion coefficient 
in Sec. IV. The possibility of anomalous diffusion in sin- 
gular topological disorder is examined by renormalization 
group theory and computer simulation in Sec. V. A dis- 
cussion of the results, and their relation to the previous 
literature, is given in Sec. VI. A discussion of the effects 
of torsion, which exists solely within the cores of defects, 
is given in Sec. VII. We conclude in Sec. VIII. 



II. THE SURFACE DIFFUSION MODEL 

We consider a particle hopping on the surface of a crys- 
tal. The particle hops only between nearest-neighbor lat- 
tice sites, and the rate of hopping is constant. In particu- 
lar, since surface diffusion is usually an activated process, 
the rate to hop between neighboring sites is assumed to 
be independent of the distance between sites. Disorder 
in the spatial arrangement of the surface lattice sites in- 
directly affects the diffusion dynamics through modifica- 
tion of the hopping events. 

We derive the Fokker-Planck, or diffusion, equation for 
the surface species by two independent methods. In the 
first method, the hopping dynamics is derived from a 
physically-motivated consideration of the master equa- 
tion for the process. In the second method, the result is 
derived in an efficient fashion by considering a change of 
variables in the field theory for the dynamics. 

The particle is considered to hop on an irregular grid 
of lattice sites. The probability for particles to be on 
a given site, P(r), decreases with time due to hopping 
of particles off the site and increases with time due to 
hopping of neighboring particles onto the site (see Figure 
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FIG. 1: A lattice site, r, on the distorted crystal and the four 
nearest neighbors are shown schematically. Also shown is the 
distorted unit cell of the central lattice site. 
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where h is the lattice spacing, At is the small time in- 
crement, and Dq is the diffusion coefficient. The volume 
of each, possibly distorted, unit cell is given by V^r). 
Equation is exact and leads in the continuum limit to 
the general expression for diffusion in curved space 
Although the crystal may be distorted, a regular crystal 
lattice can always be defined locally in terms of lattice 
coordinates <x. In the u space, the particle hops either 
up, down, left, or right. The correspondence is given by 
17 «-► <x = (7i,0), r n <-► er = {-h,0), r nl <-► a = (0,h), 
and rjv *-* & = (0,—h). The positions of the neighbor- 
ing sites are defined such that a hop in the appropriate 
direction leads to r. For example, 



r / — h 



dr 

d<7\ 



h 2 d 2 r 



2 da 2 



0(h 3 



(2) 



Note that the r coordinates are considered to be a fixed 
function of the <x coordinates: r = r(er). This mapping is 
independent of time, as the defects that generate the non- 
trivial mapping will be quenched in the two-dimensional 
crystal. Inverting eq. © for 17 gives 



. dv 

vi = v + h — - 
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With these expressions for the four neighboring sites, eq. 
JTJ to O(h) becomes 
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where the summation convention has been used. Equa- 
tion J2J is exact and leads in the continuum limit to the 



general expression for diffusion in curved space 9]. The 
notation 



9'- 



dr.i drj 
da a da a 
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will be used. The shorthand di = d/dri will also used. 
Equation is a relation for the probability distribution 
in r space. The relation for the probability distribution 
in or space requires a Jacobian: 



P(r,t) =G{eT,t)\detda Q /dri 
The Jacobian is given by 



| det da a /dn \ = </g{r) 



da x da v da x da,, 



dr„ dr„ 



dr y dr x 



(6) 
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As noted, the Jacobian is ^J~g = (det^y) 1 ' 2 = 
l/(det g n y/ 2 , where the inverse of the matrix g l i is given 
by 



da a da a 
drj dr4 



(8) 



The volume of each unit cell is given by V(r) = h 2 / -Jg{r). 
By detailed balance, since the rates to hop forward and 
back between any two sites on the crystal are the same, 
the long-time average number of particles per site must 
be equal at all sites: lim^oo G(cr, t) = const. This im- 
plies Iim t ^oo P(r, t) oc y/g{r). 
Since 
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eq. J3J becomes 



Finally, given that 
1 
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(10) 
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(see below) 



(12) 



the final expression of the diffusion equation becomes 



dC 1 



(13) 



This equation applies everywhere except within the cores 
of topological defects, because it has been assumed in the 
second to last line of eq. (|12l) that the differentiations 
commute [ijj. Equation ((13^1 is nothing more than the 
usual diffusion equation in curved space, with the familiar 
Laplace-Beltrami operator replacing the Laplacian of 
flat space. The mean-square-displacement is given by 

(r 2 (t)) = f dr|r| 2 P(r,i) (14) 



dr^(r)|r| 2 G(r,t) . (15) 



Equation (|13(l differs from the most general expression 
for diffusion in curved space by a term related to the 
torsion |9j]. We reevaluate the term d In det gj dri ■ 



1 d In det g 1 fc • d 



dn 



9jk 



2" dr. 
1 drk drj d f dap dap 



2 da a da a dn \ drj drk 
drk drj dap d dap 
da a da a drj dri drk 
drk d d 
da a dri drk 

d da a 
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+ d ^__^_\ a f 16 s 

da a \ dn dr k dr k dn J 

Defining the torsion as 2T^ k = (drk/da a )(d 2 jdridrk — 
d 2 /drkdri)a a , note that 

y 1 d In det g _ dn drj d dap ij k 
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Combining eqs. ((TUj) . (jTT|) . and (|T7jl . we find that the 
exact expression for the diffusion equation is 

dC 1 

= —di^D^djG) - 2D Q g«TidjG . (18) 



Equation Ijl8|l is equal to the general expression for dif- 
fusion in curved space The difference between the 
exact answer, eq. |JT5J, and that assuming that the or- 
der of differentiation commutes, eq. Ijl3|l . is given by the 
torsion term. The torsion is an explicit measure of the 
non-commutativity of differentiation and is, therefore, a 
measure of the defect density 01- The diffusion equa- 
tion does not apply within the cores of defects, where 
the metric tensor is undefined, and the only place where 
the torsion is non-zero. The effects of the torsion should 
probably be studied with a detailed model rather than 
with the long-wavelength, continuum theory of the dif- 
fusion equation. For this reason, we exclude this torsion 
term (although see section VII below). The long range, 
external to defect core, effects of the topological defects 
are, of course, included in eq. ifH^ through the metric 
tensor <j> y and *Jg~. Refs. [a, 0] included the torsion term 
explicitly, and a series of approximations allowed the gen- 
eration of non-physical dynamics. 

Equation (|13f) can, alternatively, be derived by consid- 
eration of the field-theoretic representation of the diffu- 
sion operator [l2l Il3[ . In this representation, the Green 
function is given by an average over a field: 



G(r,i) = («M)) , 



(19) 



where the average is taken with respect to the weight 
exp(— S). The particle hopping occurs in a space with- 
out regard to the distortion of the crystal, as the rate of 
hopping is independent of the distance between lattice 
sites. The action for such normal diffusion is given by 

"5 = J dtj da- japt + 6(t)]a - A>a|^rj 

+ / d<xn o (<T)a(r,0) ," (20) 



where no is the initial density profile, and details of the 
replica indices used to accommodate averaging over dis- 
order have been suppressed 0, . This action is en- 
closed in quotations since the er space is not well-defined 
in the presence of topological defects. That the diffu- 
sion is normal in cr space, however, does make it clear 
that the limiting distribution should be limt_»oo G(er, t) = 
(const). From eq. ©, then, the limiting distribution in 
r space is given by limt_>oo P( r , t) = (const)^/p(r) = 
y/g(r)/ Jdr'y/g(r'). While this result may be surpris- 
ing, note that the defects which distort the geometry 
must affect the limiting distribution, unlike the typi- 
cal case in differential geometry where the observables 
are described by a theory independent of the coordinate 
system. This explicit result for the limiting distribution 
agrees with the prediction from the simple detailed bal- 
ance argument given above. Note that J dY^fg(r) is a 
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constant for a given realization of the quenched disorder. 
The long-time normalization factor for the probability is 
fixed to be the inverse of this integral by the initial con- 
dition P(r, 0) = S(r). Equation (|13fl for the dynamics 
conserves J dr^/g(r)G(r, t), hence, the probability distri- 
bution, P(r,t) — yfg(r)G(r,t), remains normalized to 
unity for all times t > 0. After change of variables from 
<t to r, again making the assumption of being outside 
defect cores so that differentiation commutes, the action 
becomes 



S = J dt I dv{^ga[d t + 5{t)]a-ad t [^gD^d 3 a]} 



+ / dr^/gno(r)a(r, 0) . 



(21) 



Finally, integrating out the a field, using eq. 1|19|1 . and 
noting that for the Green function rio(cr) = (5(er), the 
Fokker-Planck equation is 



(22) 



with G(r,0) = S(r)/ \fg(r). The field-theoretic result, 
eq. ill' 21) . is the same as that derived by more physically- 
motivated means, eq. (I13f) . 



III. THE MODEL OF TOPOLOGICAL 
DISORDER 

The topological defects modify the diffusive motion of 
the particle by affecting the g^ in the Fokker-Planck 
equation. Once g^ is determined, eqs. (fT3|) and i(T5|) 
provide the means to calculate the transport properties. 
It is conventional in continuum elasticity theory to relate 
the spatial coordinates to the lattice coordinates by 



r(cr) 



u(r) 



(23) 



where the displacement field u is written in terms of the r 
variables that remain well-defined even in the presence of 
topological defects. The <x space, on the other hand, does 
not remain well-defined, since the effect of disclinations 
is to add or remove wedges of lattice sites from er space, 
and the effect of dislocations is to add or remove half- 
lines of lattice sites from er space. For a dislocation at 
the origin with Burgers vector b, the displacement fields 
are given by [l6l | 



2nu 



disloc 



(/! + A) e k ibmr k _ x r y 

+ bi tan — 



(2/u + X) r 2 

eubi in — 

2/i + A h 



(24) 



where \x and A are the two-dimensional Lame coefficients, 
e n = e 22 = 0, and e 12 — —£21 = 1- Similarly, for a 
disclination of strength s at the origin, the displacement 



fields are given by 
27raf sclin = 



(M + A) 



2(2^ + A) 



se lk r k tan" 



/J 



2^ + A 



sri \n{r/h) 



(25) 



Equation l|25|) differs from the simplified distortion field 
used in by the inclusion of the strain field representing 
the local lattice contraction and expansion. These are the 
terms in eq. (I25|l that depend on the Lame coefficients. 
Since linear elasticity theory is used, the dislocation field 
is given by the dipole limit of two superimposed discli- 
nation fields: 



disloc 



(bi/s)e 3 id 3 u. 



disclin 



const . 



(26) 



The derivatives of the displacement fields are required 
to evaluate g^ from eq. ©. The dislocation fields are 
preferable for this calculation, as they lead to well-defined 
Fourier transforms: 



disloc 



A 2k 2 x k y 



fc 2 



+1 



2/i + A fc 4 

/i + A k x {-k 2 x + kl) 



6x(k) 



OyU X 



disloc 



2^ + A 

^ + A 2k x k, 



fc 4 



2/j. + A k 2 



by(k) 



k 2 



+1 



disloc 



2/i + A k 4 

/i + A ky{-kl + kl) 
2/i + A /fc 4 

/i + A k x {-k 2 x + k 2 y ) 



&x(k) 



2/i 
/i 



Xk 2 
k 



+ 1 



2/i + A 
/i + 



fc 4 



2/i + A k 2 



6 v (k) 
fex(k) 



ByUy 



disloc 



2/i + A k 4 h 
/i + A k y (-k 2 x + k 2 y ) 



by(k) 



// 



k,. 



2/i + A 



k 4 



2/i + A k 4 



k x 
k 2 



2/i + A k 2 

by(k) . 



6x(k) 
(27) 



Note that the x and y derivatives of the strain fields are 
not simply related by the ratio k x /k y , due to the pres- 
ence of the defects. The linearity of elasticity theory has 
been used to accommodate a density field of defects with 
Burgers vectors given by b(r). 

The dislocations are assumed to be distributed ran- 
domly in the material with correlation function 



= ^(27r) 2 ,5(k 1 + k 2 )x(ki + k 2 ) , (28) 



where 



X(k) =jk n exp(-f3k 2 ) 



(29) 
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Physically, we expect this model of dislocations should 
generate identical dynamics to one in which disclinations 
are randomly distributed with correlation function 

(s(k 1 )J(k 2 )> = (2 7 r) 2 ^(k 1 +k 2 )|ki +k 2 | 2 x(k! +k 2 ) . 

(30) 

This physical expectation is a mathematical consequence 
of eq. PI . 

With these results in hand, we are now in a position to 
calculate the action for the field theoretic description of 
the Green function. The terms in eq. (|21|) are expressed 
to linear and quadratic order in it dlsloc 1 and then an aver- 
age over the random distribution of dislocations is taken. 
In fact, since eq. Ijl4(l is preferable to eq. (|15[) . the theory 
is written in terms of the fields a, c, where c = y/ga, and 
P = (c). The action is 

S = Jo dt I dT ^ dt ~ + SDl1 ^ 2 + S ^ C 

+ J dr8(r)a(r, 0) + Si , (31) 

where 

= -2D 2 [ dt x dt 2 



'kik2k 3 k4 

x(2 7 r) 2 <5(k 1 +k 2 +k 3 +k 4 ) 
xa(ki, ti)c(k 2 , ti)a(k 3 , i 2 )c(k 4 , t 2 ) 



fi + X ki -k 2 (fc 2 + fc^ 



2[i + \ 1 2(i + A 



Iki+^P 



H 2 fj, + A k 3 • k 4 (fcf + fcf) + 2fc 2 fcf 



2^ + A J 2^ + A 
X(ki + k 2 ) 
Ikr+k,! 2 ' 



|k 3 +k 4 | i 



(32) 



where the notation J k stands for J" d 2 k/(27r) 2 . The term 
resulting from a non-zero average of (u dlsloc ) 2 is 



SD 



ii 



2tt 



13A 2 + 16/xA + 25/i 2 



9(2^ + A) 2 
dk fc"- 1 exp(-/3fc 2 ) 



7 A) T(n/2) 
2tt 2/W 2 



13A 2 + 16/xA + 25^ 2 



9(2 M + A) 2 



(33) 



Exactly the same theory is generated if the correlation 
function eq. I|3U|I is used with the disclination displace- 
ments given by eq. (|25|l . 



IV. TOPOLOGICAL DISORDER REDUCES 
THE DIFFUSION CONSTANT 

For the model with n > 0, the topological disor- 
der reduces the diffusion coefficient by a finite amount. 



The finite contribution of SDn is explicit in eq. I|33|l . 
Moreover, standard power counting arguments [lTj show 
that non-perturbative, renormalization effects can be ex- 
pected from eq. (|32|) only for n < 0. From perturbation 
theory on eq. (|32() for n > 0, the contribution to the 
diffusion coefficient is found to be 



5Dr 



7 A) r(n/2) 
' 2tt 2/W 2 



4/i 2 + 2(/i + A) 2 



(2^ + A) 2 



(34) 



The total contribution to the diffusion coefficient is, 
therefore, 



SD = - 



7 A) r(n/2) 
2tt 2/W 2 



29/r 2 + 20/zA + 5A 2 



9(2^ + A) 2 



(35) 



To demonstrate the behavior of this model, we perform 
numerical simulations. The dislocation density fields are 
constructed with correlation function eq. I|28|) for n = 2 
using the method of ref. [l8|. Equation Q27[l and an in- 
verse fast Fourier transform are used to calculate the dis- 
placement fields in real space. The matrix g lJ is calcu- 
lated as the inverse of the matrix gij given by eq. JSJl, 
and the relation cr = r u(r) is used. 

The Fokker-Planck equation, eq. (|13|l . can be consid- 
ered to result from many small hops, the net effect of 
which is Gaussian, diffusive motion. So that a hopping 
process on a lattice reproduces this differential equation, 
the average and mean-square displacements must be cor- 
rect at each lattice site. Interestingly, this differential 
equation can be evaluated by Monte Carlo methods on a 
perfect, square lattice, even though the differential equa- 
tion itself describes the motion of a particle in a distorted 
geometry. To first order in the time step, the mean dis- 
placement is given by 

(n(At)) = J dr^g ri G(r,At) 

= J dt J dv^g nd t G(T,t) 



(36) 



dr^g n G(r,0) 
D At. 



V9 



-9AV99 1 



where eq. Ijl3|l and integration by parts twice has been 
used in the last step. Similarly, to first order, the mean- 
square displacement is given by 

{n{At) rj (At)) = J dr^gr irj G(r,At) 

= J dt J dr^/g nrjd t G{T,t) 

+ J dr^ r^ GM) 
= 2D Atg i i , (37) 

where eq. (|13fl and integration by parts twice has again 
been used in the last step. 
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Eight hopping rates are defined, consistent with the 
specifications of eqs. <|5oTJ!T7J) . So that the non-diagonal 
terms of g lJ are properly reproduced, both nearest- and 
next-nearest-neighbor hops are required. The rate for 
each hopping event is 



T,(r 



Ar) 



h 2 



g(r + Ar) 



-l 1/4 



S(r) 



x-[/(r)+/(r + Ar)] 



(38) 



The function / is given by 



/ 
/ 
/ 
/ 



= g 11 - e, for Ar = (±h, 0) 
= g 22 - e, for Ar = (0, ±h) 



(g 12 

(-9 



e)/2, for Ar = (±h,±h) 
12 1 e)/2, for Ar = (±h,^h) 



(39) 



with e = |(? 12 |. The transition rates in eq. (|38|) explic- 
itly satisfy detailed balance for the equilibrium distribu- 
tion lim t _ >00 P(r, t) = (const) ^/gYr). These rates give the 
correct average and mean-square displacements to 0(h), 
eqs. (|36H37|) . when At = l/(Ei T i)- Th ese results imply 
that the Monte Carlo procedure evaluates the differen- 
tial eq. and so eq. (f 141) can be used to calculate the 
mean-square-displacements. The procedure of ref. [18l | is 
used to perform the simulation of this random process, 
where the particle is moved to one of the neighboring 
eight sites with probability AtTi, and time is incremented 
by dt = —At ln(x), where a; is a uniform random number, 
< x < 1. 

The results of the numerical simulations are shown in 
Fig. [21 The calculations were performed for the case /i = 
A, n = 2, h = 1, and (3 = 4. The simulations were 
done on 4096 x 4096 lattices for a total of 500000 steps 
and averaged over 100000 particles. The strength of the 
disorder was varied between < 7 < 1.25. For larger 
values of 7, the transition rates specified by eqs. (|3^1 - l3^|) 
became negative at some of the lattice sites. Also shown 
is a fit to the functional form SD/Dq = 1 — ax. The fit to 
the simulation data of a = 0.01285±0.0051 is in excellent 
agreement with the theoretical of result a = 0.01326 from 
eq. lEHJ- 



V. POSSIBLE ANOMALOUS DIFFUSION 

The case n = is interesting, as perturbation theory 
for the diffusion coefficient formally diverges. While this 
theory has the same upper critical dimension, d c = 2, 
as the problem of diffusion of an ion in the electro- 
static field of random, quenched charges [l|, the inter- 
action term, eq. I|32l) . is quite different. In comparison 
to the analogous term for diffusion in the random po- 
tential (e.g. term S3 of ref. 0] with x vv (k) — x(k)/k 2 ), 
the term proportional to \x is new, as are the factors 
-k\ + 2[k\k\ - (ki • k 2 ) 2 ]/|ki + k 2 | 2 in the term pro- 
portional to /1 + A. Indeed, as we will see, the present 
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FIG. 2: Shown are simulation results for the reduction in the 
diffusion coefficient for the case n = 2, fj, = A, and /3 = 4. 
The error bars are roughly ±0.01. The best linear fit to the 
simulation data is shown (solid line). The simulation data 
are compared to perturbation theory (dashed line), eq. 1351 
5D/D = -7/(6vr/3). 



interaction term is more difficult to analyze than is the 
analogous one from diffusion in a random potential. For- 
mally, the case of n < leads to large distortions of the 
lattice for arbitrarily small 7, which implies that the as- 
sumption of linear elasticity used to calculate the strain 
fields breaks down. We can, however, treat the dynamical 
behavior implied by eqs. I|31H32(I as an interesting math- 
ematical question. A technical detail is that we supple- 
ment the correlation function eq. (|29|l with the condition 
x(0) = so that the displacement fields of eq. (|2*7Jl are 
well-defined for k = 0. This suppresses macroscopic size 
fluctuations of the sample. 

Before applying renormalization group theory, the 
terms in the field theory must be known. The quartic 
interaction term, eq. I|32|l . is known. The contribution to 
the propagator, eq. (|33|l . while explicit, leads to a formal 
divergence of the short-time diffusion coefficient. Nu- 
merical simulations show that the local diffusivity tensor, 
Dog IJ , can be large but is never vanishingly small. The 
locations of large local diffusivity, moreover, are isolated. 
The apparent divergence of 5Dn is, thus, simply the re- 
sult of particles rapidly hopping away from a few isolated 
locations. These physical considerations suggest that the 
divergence of SDjj is washed out by spatial averaging 
and is not important for the long-time dynamics. We 
can, therefore, assume a finite local diffusivity. Numeri- 
cal simulations of the dynamics, to be described below, 
bear out this assumption of a finite short-time diffusivity. 
Indeed, a finite short-time diffusivity is assured for finite 
lattice sizes by the elimination of the mode. The 
anomalous dynamics, then, is observed on finite lattices 
for time scales that are less than the characteristic time 
it takes to travel across the lattice. 

We apply renormalization group theory to the action 
(|31H32|) to take into account the effects of nonzero 7. To 
one-loop order, self-energy and vertex diagrams are sum- 
marized in Figs. |21 and 2| The flow equations are inte- 
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FIG. 3: a) Diagram representing the propagator. The arrow 
points in the direction of increasing time, and double lines 
represent the bar fields, b) Disorder vertex 7. 







FIG. 4: One-loop diagrams: a) self-energy diagrams con- 
tributing to D. b,c,d) vertex diagrams contributing to 7. 
Diagrams (c) and (d) cancel. 



grated to a time small enough so that perturbation theory 
applies. In this regime, matching theory is used to deter- 
mine the constants of integration for the flow equations. 
Momenta in the range A/6 < k < A are integrated over, 
and the fields are rescaled by a'(6k, b~ z t) — a(k,t)/a 
and c'(Mc, b~ z t) — c(k, t)/a. The relations a = 1, a — b 2 
are used to achieve a fixed point and to keep the time 
derivative in S constant. The flow parameter is defined 
by I — ln6. We determine the dynamical exponent, z, by 
requiring that the diffusion coefficient remain unchanged. 
Defining 7l 2 = 1 ^ 2 /(2^+X) 2 , 7 2 = 7A1 ( M +A)/[2(2 M +A) 2 ], 
and 73 = 7(/i + A) 2 /[4(2/i + A) 2 ], the contributions to the 
parameters from the one-loop diagrams of fig. Q] arc 



dln 7l 2 

dl 
din 7I 

dl 

dl 



2(z-2)--(7 2 + 27 2 ) 

2(z- 2)--(7 2 + 37 3 2 - 27i73) 
2(^- 2)--(73 2 - 2 7l 73) . 

7T 



(40) 



From the requirement that the diffusion coefficient re- 
main fixed, the dynamical exponent is 



= 2 + -(7 2 + 2 7 2 ) . 



(41) 



Using eq. I|41|) in cq. 1)40(1 . the flow equations become 



din 7 2 

dl 
dhi7 2 

dl 
dhi7f 

dl 



(71 +73) 2 



-(71+73)' 

7T 



(42) 



As expected, the flow equations show that there are only 
two independent parameters, 71 and 73, resulting from 
renormalization of the two Lame coefficients. In other 
words, the relation 7! (0 = 7i (073(0 I s maintained under 
the renormalization. 

Unexpectedly, however, the flow equations show that 
the 72(0 an( i 73(0 are growing. Indeed, these one-loop 
flow equations predict 73 (Z) flows to infinity at a finite 
time corresponding to / = [7r/(27°) 2 ]{ln[(7] ) + 73O/73] — 
7i/(7i +73)}- The divergence of this parameter implies 
that higher order terms must be kept in the flow equation 
to derive a controlled result. It may also be the case 
that terms higher order in u dlsloc must be kept in the 
expansion of the action l|21|l. 

If the renormalization of the parameters is assumed to 
be controlled by higher-loop corrections and small, the 
dynamical exponent can be used to determine the scaling 
exponent for the mean-square displacement at long times: 



(r 2 {t)) - (const)* 



1-5 



The renormalized time flows as 



t(l*)=te~fo z{l)dl 



to 



(43) 



(44) 



where the flow equations are stopped at I* so that *o ~ 
h 2 / (AD ). The renormalized mean-square displacement 
flows as 



(r 2 (t))=e 2l (r 2 [t(l),l}) . 
Finally, at the matching 

(r 2 [t(ni*])=4Dt(l*) , 



(45) 



(46) 



since the time is short enough so that the disorder does 
not significantly affect the motion of the particle. In 
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other words, it is assumed that at short times the diffu- 
sion coefficient remains finite, despite the formal appear- 
ance of 5Dij in eq. . Putting these matching results 
together with the dynamical exponent, the mean-square 
displacement is found to scale at long times as 

(r 2 (i)) ~ (conBt)t 1 /[ 1 +(^+»rg)/*] . (47) 

To test whether anomalous scaling occurs in the 
full non-linear model, we perform numerical simula- 
tions. The transition rates from eqs. (|38H39[1 cannot 
be used, as they are negative even for small values of 
7. We, therefore, develop a new strategy based upon 
the idea that diffusion locally follows a Gaussian prob- 
ability distribution with mean and variance specified by 
eqs. H36I - I37JI . The time increment At is chosen so that 
max(|(5rj)|, \{8ri5rj)\) is on the order of unity. This is 
done by choosing 1/At to be the maximum of the abso- 
lute values of the two average displacements in eq. jSEjl 
and the two eigenvalues of the matrix D$g l i . Defining the 
matrix w = (2DoAtg i: ') 1 ^ 2 , the random displacements of 
the diffusing particle are given by the relations 

Ay = ^tyy/gg*) + w n Zl + w 22 z 2 t (48) 

where Z\ and z 2 are independent, Gaussian random vari- 
ables with zero mean and unit variance. This approach 
reproduces the Fokker-Planck equation (0J in the limit of 
a small lattice spacing and time increment. For a finite 
lattice spacing, the diffusion coefficient, and possibly the 
scaling exponent 5, contain discretization errors. 

The results of numerical simulations with this scheme 
are shown in Fig. [5] The calculations were performed for 
the case /u = A, n = 0, h = 1, and j3 = 4. The simulations 
were done on 4096 x 4096 lattices for a total of 1000000 
steps and averaged over f 00000 particles. Also shown is 
a fit to the functional form of eq. I|43(l . The simulation 
results are approximately fit by (z — 2)/j — 1.02 ± 0.28. 
If it is assumed that none of the parameters flow, the 
scaling exponent is given by (z~2)/j — 2/(3ir) 0.2122. 
By comparison, the simulation results suggest that there 
is substantial positive renormalization of the parameters, 
as is suggested by eq. (@2J>- 

VI. DISCUSSION 

A Fokker-Planck equation for diffusion on the surface 
of a crystal with topological defects, eq. has been 

derived by two independent methods. As expected, the 
usual diffusion equation in curved space is derived. An 
additional assumption of ^fg = 1 of previous, approxi- 
mate treatments |6j has also been removed in the present 
calculation through the use of the exact strain field eq. 
1)25(1. The theory of random dislocations is shown to be 




Y 



FIG. 5: Shown are simulation results for the scaling exponent 
for the case n = 0, fj, = \, and f3 — 4. The error bars are given 
roughly by the scatter in the data. The best linear fit to the 
simulation data is shown (solid line). Also shown (dashed 
line) is the prediction assuming that the other parameters do 
not flow, z = 2 + 27/(3-71-). 



equivalent to a theory of random disclinations, where a 
simple factor of k 2 relates the correlation functions of the 
two models of disorder, eqs. I|28(l and (|30|l . 

The field theory for disorder, eqs. (|3 1 H32fl . is explicitly 
shown to be distinct from that for diffusion of an ion in 
a random electrostatic potential field. One consequence 
of this difference is that the renormalization group flow 
equations are more involved to analyze, with one-loop 
results unable to render a controlled prediction. 

Topological disorder slows down a diffusing particle, as 
shown by eq. (|35|) . This reduced transport should be ob- 
servable on the surfaces of crystals with quenched discli- 
nation or dislocation defects. While the effect is subtle, 
it would be an interesting one to observe experimentally. 
The present computer simulation results suggest such ob- 
servations should be feasible. 

For singular disorder, n < in two dimensions, the 
model of topological disorder leads to subdiffusive motion 
of the particle. Of course, for such singular disorder, the 
assumption of linear elasticity breaks down. Moreover, 
the energy of a distribution of topological defects with 
net dipole moment becomes super-extensive due to large 
strain fields at the edges of the two-dimensional crystal 
[To|. Nonetheless, the suggestion that subdiffusion is the 
mathematical result of motion in the, possibly approxi- 
mate, random displacement fields of linear elasticity the- 
ory is interesting. Renormalization group arguments are 
suggestive of such subdiffusion, although one-loop results 
are unable to capture the exponents quantitatively. 

Numerical simulations accurate to all orders in the dis- 
placement fields suggest that the motion is, indeed, sub- 
diffusive. These numerical simulations suggest that there 
is significant renormalization of the disorder strength pa- 
rameter, in contrast to the case of diffusion in random 
potential fields 0. Interestingly, the renormalization of 
z appears less significant for smaller values of 7, although 
this may be because the crossover time for renormaliza- 



9 



tion is large for small 7 and longer than the observed 
simulation time. These simulations suggest a power law 
behavior of the mean square displacement, although lo- 
calization at exceptionally long times cannot be ruled 
out, in principle. 



VII. WHAT IF TORSION IS INCLUDED? 

We here comment on the impact of the torsion term 
within the continuum theory of the diffusion equation. 
The torsion term is evaluated as 



\da a )\dndr k dr h i)v, 



(49) 



Expanding eq. (|18() to linear order in it dlsloc , we find that 
the interaction term, previously eq. I(32|l . becomes 




FIG. 6: Shown are simulation results for the scaling exponent 
for the case n = 0, \i = A, and /3 = 4 when torsion is included. 
The error bars are given roughly by the scatter in the data. 
The best linear fit to the simulation data is shown (solid line). 
Also shown (dashed line) is the prediction assuming that the 
other parameters do not flow, z = 2 + 47/ (3tv). 



Si = 



D 2 f 

— / dhdh 

* J Jkik 2 k 3 k 4 

x(2^) 2 5(k 1 +k 2 +k 3 + k 4 ) 
xa(ki, ii)c(k 2 , £i)»(k 3 , t 2 )c(k 4 ,t 2 ) 



20 



fcf 



47: 



ki • k 2 (k 2 + fcf) + 2k\k\ 
|ki+k 2 | 4 



x(k x 



271 



^3 



k 4 P 



-4 7 . 



k 3 • k 4 (fc 



x(k 3 + k 4 ) + v ^k 4 
xx(ki +k 2 ) . 



(50) 



Exactly the same theory is generated if the correlation 
function eq. 1)30(1 is used with the disclination displace- 
ments given by eq. ((25(1 . The inclusion of the torsion term 
has generated the additional terms proportional to y/jk 2 
and 1/7 A: 4 . Applying perturbation theory to Si, we find 
that a mass term is generated, 5m — — 2Do7iy / 7~. This 
term is exactly canceled by a mass term arising from the 



average of terms proportional to 



.disloc\2 



) , which must 



be the case since the master equation Q conserves prob- 
ability. No contribution to the diffusivity is generated by 



the average average of terms proportional to (u dlsloc ) 2 . 
From the average of Si, we find an additional negative 
contribution to the diffusivity: 



SDin = - 



7A> r(n/2) 
2tt 2/W 2 



2(^ + A) 
(2^ + A) 



(51) 



Within the approximation of the continuum diffusion 
equation, then, the torsion term generates an additional 
contribution to the effective diffusivity when n > 0. Note 
that this contribution is a result of correlated drift terms 



that exist solely within the cores of the defects. There is 
no reason to expect that this contribution is universal or 
even well-described by continuum theory. 

For the mathematically interesting case of n < we 
follow our previous numerical strategy. The random dis- 
placements of the diffusing particle are altered from eq. 
6HJ) to 



Ax = 



D At 



V9 

(52) 

To make use of this formula, we need an expression for 
dri/da a that occurs in Tj. This is found as dri/da a = 
where A a i = da a /dri = 5 a i — du a /dri. In evalu- 
ating Tj, we use the first line of eq. ((49(1 . The results of 
numerical simulations with this scheme are shown in Fig. 
El The calculations were performed for the case [i = A, 
n = 0, h = 1, and /3 = 4. The simulations were done 
on 4096 x 4096 lattices for a total of 1000000 steps and 
averaged over 100000 particles. Also shown is a fit to the 
functional form of eq. (|43|) . The simulation results are 
approximately fit by (z -2)/^= 0.49 ± 0.04. 

There appears to be relatively little if any renormal- 
ization of z away from the bare value. A power law be- 
havior of the long-time mean square displacement in the 
presence of torsion is observed, although localization at 
exceptionally long times still cannot be ruled out. 

It is clear that within the continuum assumption of the 
diffusion equation, the torsion term affects the dynamics. 
The contribution to the diffusion coefficient is explicit in 
eq. 1(51(1 for the case n > 0. For n = 0, the results 
shown in Fig. differ from those without torsion in Fig. 
03 Note that the results with torsion, as those without 
torsion, differ substantially from the approximate results 
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of |a, l3 1 noticeably through their dependence on the two 
Lame coefficients. 



VIII. CONCLUSION 

We have given a treatment of the effect of topolog- 
ical disorder on transport properties. Within the lat- 
tice reconstruction predicted by linear elasticity theory, 
topological disorder is manifestly different from charged, 
potential-type disorder. The net effect of the defects, 
through local lattice expansion and contraction and 



global topological rearrangement of lattice connectivity, 
is an overall reduction of the transport. Interestingly, 
randomly placed dislocations, or randomly placed discli- 
nations with no net disclinicity, lead to anomalous sub- 
diffusive behavior when the displacement fields of linear 
elasticity are used. 
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